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^3" (57) Abstract: A signal processing utility is disclosed involving time-to-frequency domain transforms for applications including 
O medical diagnostic signal processing. Such transforms can be used to define a continuous spectral density function or other spectral 
^ density function including non-zero values at irregularly spaced frequency intervals. The invention thereby enables more accurate 

representation of certain real spectra and reduced spectral broadening. The utility also accounts for digitization errors associated with 
Q analog-to-digital conversion. The invention has particular advantages with respect to medical contexts where the received signal has 

a changing spectral content as a result of interaction of an interrogating signal with moving physiological material such as blood 
^ flowing through an artery. 
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DIAGNOSTIC SIGNAL PROCESSING METHOD AND SYSTEM 

FIELD OF THE INVENTION 

The present invention relates generally to signal processing applications and, in particular, to 
signal processing applications involving time-to-frequency domain transforms where it is useful to 
define a spectrum including frequency information for frequencies other than a fundamental 
5 frequency and harmonies thereof, or where it is useful to generate a spectrum based on a portion of 
the time-based signal representing less than a period of a signal component. The invention has 
particular application to methods wherein a received signal is decomposed or otherwise processed 
in a manner that provides information regarding properties of a material modulating or otherwise 
producing the received signal. The received signal may emanate from tissue interacting with an 
10 interrogating signal that was introduced into an organism and was modified from the interrogating 
signal by one or more influences of the organism's material, e.g., organs, tissue, fluids or other 
structure. 

BACKGROUND OF THE INVENTION 

15 Noninvasive diagnoses using interrogating signals are frequently employed to ascertain the 

condition of tissue and functional parameters of organs. Such interrogating signals may use, for 
example, electromagnetic or sonic, including ultrasonic, energy. One example is transmitting 
ultrasonic waves into an organism, such as a human being, and receiving a return signal that is a 
modified version of the transmitted interrogating signal, with the modifications caused by the 

20 properties of one or more tissue types. The received signal is measured and processed to evaluate 
the condition or properties of the tissue(s) or other material(s) with which the interrogating signal 
interacts. In particular, interrogating ultrasonic signals are often transmitted into an organism and the 
received signal is due to backscattering or echos of the interrogating signal by one or more tissue 
types. Even more particularly, an ultrasonic interrogating signal is often transmitted into a region of 

25 the organism where blood flow is expected and the received signal is a modified version of the 

interrogating signal wherein such modifications are due in whole or part to frequency shifts caused 
by Doppler effects caused by movement of tissues, including movement of blood. In this case it is 
desirable to measure and process the received signal in a manner such that the amount of Doppler 
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frequency shift that occurs can be ascertained in order to estimate the various velocities of the 
moving tissue, particularly the velocities of the flowing blood. 

Ultrasonic Doppler techniques typically entail using a hand-held probe that contains one or 
more transducers that are used to transmit the interrogating signal and to produce a response, 
5 typically in the form of a time- varying voltage, to the received signal. The interrogating ultrasonic 
signal typically has a frequency in the radio frequency (RF) range of about 0.5 million Hertz (0.5 
MHz) to about 20 MHz, although other frequencies may be used. The interrogating signal may be 
sent continuously, in which case the method employed is termed "continuous wave" (CW) or the 
interrogating signal may be sent intermittently, usually on a periodic basis, in which case the method 

10 employed is termed "pulsed." When pulsed Doppler is used it is possible to interrogate tissue over 
a range of specified distances from the transducers by measuring the received signal only for a 
selected duration that occurs after a predetermined time delay after the pulse is sent. If the speed of 
sound for the intervening tissue is known then the time delay converts directly into a depth of tissue 
being interrogated. Interrogating a specified tissue depth often provides clinically valuable 

15 information that cannot be obtained if CW Doppler is used. CW Doppler systems are often simpler 
than pulsed Doppler systems, although they cannot provide information that is known to come from 
tissue associated with a specific time delay (and corresponding estimated tissue depth). 

Regardless of whether the Doppler system is pulsed or continuous, the spectral content of 
the received signal is analyzed to assess how the tissue modified the interrogating signal. The 

20 spectral frequency content of the signal is a measure of the amplitudes of each of the constituent 

frequencies that, when combined, form the total signal. The spectral content is typically represented 
as either the amplitude of each of the frequencies or as the power of each of the frequencies. When 
the spectral content data are graphed the resulting graphs have a variety of names including 
spectrograms and spectral density functions. The term "spectral density function" (SD function) will 

25 be used as the general term below. When the frequency content of the interrogating signal is known 
it can be compared to the spectrum of the received signal and the differences between the spectra 
can be interpreted as arising from the velocity profiles of the tissues through which the signals 
passed. The velocity profiles can be back-calculated from spectral differences using established 
methods that are well known. 
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To date, a variety of methods have been used to estimate the spectral frequency content of 
the received signals. The main spectral analysis methods are traditionally classified as being either 
nonparametric or parametric methods. 

The nonparametric methods are related to using Fourier series and they estimate the 
5 spectral content at equally spaced frequency points, The extant methods include Fourier and 

z-transform (including the Chirp z-transform). For brevity, all methods which employ assessing 
spectral content at equally spaced frequency points are referred to as Fourier methods. 

Fourier methods, typically in the form of fast Fourier transforms (FFT), are well known and 
decompose a signal into frequencies that are integral multiples of a fundamental frequency. The 
10 integral multiple frequencies form an orthogonal basis in the frequency domain. This information can 
be used to approximate the signal as a trigonometric polynomial, as is well known. Fourier methods 
have been widely used to analyze the received signal because such methods have been refined to 
the point where they are fast and have well defined properties. 

The other class of spectral analysis methods are parametric methods. Parametric methods 
15 use time-series analysis to estimate the parameters for a rational function model. The rational 

function can have either only poles, which is an autonegressive (AR) model, only zeros, which is a 
moving average (MA) model, or both poles and zeros, which is an autoregressive-moving average 
(ARMA) model. Included in these methods are maximum entropy methods (MEM). The number 
of model parameters used in a parametric method is termed the order of the model. The various 
20 methods that exist for calculating the values of the model parameters, such as the Burg algorithm 
and Durbin's first and second methods, are well known. An advantage of parametric methods is 
that they can produce reasonable spectral estimates using very few cycles or fractions of cycles of 
data. 

25 SUMMARY OF THE INVENTION 

It has been recognized that alternatespectral analysis methods including alternate time-to- 
frequency domain transforms may be beneficial for a variety of signal processing applications, 
including certain medical applications. In particular, in spite of the extensive use of Fourier methods, 
they suffer from significant deficiencies. Parametric methods also suffer significant deficiencies. 
30 With regard to Fourier methods, the most significant deficiency is that these methods produce 
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spectral estimates for only a small number of discrete frequencies, even when it is expected or 
known that the real spectrum covers a continuous range of frequencies or that the spectrum may 
contain frequencies different from those that form the basis that is the output of the Fourier 
transform For example, the signal may contain frequencies that are not an integral multiple of the 
fundamental frequency. When the basis is not correct, either because the underlying spectrum is 
continuous or because it contains frequencies other than those in the basis of the output of the 
Fourier transform, significant errors are produced. These errors manifest themselves as a spreading 
of the spectrum in which amplitude (or power) is spread out over frequencies that are not actually in 
the signal. Fourier methods approximate the actual frequency content of the signal by spreading 
amplitudes to frequencies away from the basis frequencies. Such spreading is a well-known 
problem with Fourier methods and is called spectral broadening. Spectral broadening cannot be 
eliminated, although various methods have been employed to reduce it such as those that rely on 
windows. These windowing methods are not relevant to the current invention other than their well 
known existence and extensive use confirms that spectral broadening is a significant problem caused 
by the discrete approximations that are inherent in Fourier methods. Spectral broadening will 
always occur when discrete frequencies are used to estimate the spectral content of signals that 
have spectra that are either continuous or composed of such a multitude of closely spaced 
frequencies that they are essentially continuous, such as the signals generated by Doppler shifting of 
an interrogating signal from blood flowing in an artery. 

Another major problem of Fourier methods is that they cannot be used to examine all 
frequencies of interest. They can only include frequencies that are harmonics of the fundamental 
frequency used in the measurement. These problems, again, are well known and they limit the 
ability of Fourier methods to examine frequencies that are lower than the fundamental frequency 
used in the measurement Furthermore, this limitation leads to the frequency resolution equaling the 
fundamental frequency. For example, to estimate the amplitudes of frequencies spaced 20 Hz apart 
the fundamental frequency must be 20 Hz or lower. Therefore, if an interrogating signal has a 
frequency of 2 MHz and the need exists to discriminate the received signal with a resolution of 20 
Hz than at least 100,000 harmonics are in the basis. This large number of harmonics creates 
significant computational difficulties and is generally regarded as impractical for a variety of reasons, 
including the requirement to use many data points to perform such calculations. 
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The limitation that all basis frequencies for Fourier series are based on a fundamental 
frequency also causes other problems, such as when the received signal is not stationary (varies with 
an independent variable, such as time). This situation will exist, for example, when the received 
signal has its spectral content changing as the result of the interrogating signal interacting with blood 
5 flowing through an artery. The velocity profile changes in an artery as the blood pressure varies in 
response to the beating of the heart and the variations in the blood's velocity profile will lead to a 
nonstationaiy signal that has a time-varying SD function. If the SD function is to be measured over 
a short time interval, such as over 0.01 second, then the fundamental frequency for the Fourier 
analysis cannot provide a better frequency resolution than that which is associated with the 

10 measurement time interval. For example, a data sample with a period of 0.01 second will have a 

fundamental frequency of 100 Hz, which precludes calculating the spectral content of signals with a 
frequency resolution better than 100 Hz . 

Parametric methods also suffer deficiencies that are significant for certain applications. In 
particular, these methods minimize the sum of the squared errors between the measured data and 

15 the values estimated by the model. As described later, measured data are necessarily digitized and 
this process introduces errors that produce a wide variety of equally valid models. Minimizing least 
squared errors does not account for the possibility of such alternative solutions. 

Another shortcoming of parametric methods is their critical dependence on the order of the 
model chosen. For example, in the most common model, AR models, if too few poles are used 

20 then spectral peaks are not reliably detected. If too many poles are selected then false peaks are 
introduced. Proper model selection is nontrivial and requires a computationally intensive search 
process that involves calculating models using a range of orders and then selecting among them. 

It has also been recognized that digitization errors can significantly affect spectral analyses in 
certain applications. The input analog signal is not used directly to calculate spectral density 

25 functions. First, the measured data are digitized and given discrete values. This process quantizes 
the data and introduces errors. These digitization errors (also called discretization errors and 
quantization errors) confound spectral analysis. Each digitized measurement is an estimate of the 
actual value and has an upper and lower error limit, within which any value is equally likely to occur. 
Even though any value within the error limits is equally likely to be the true value, the digitized value 

30 is conventionally established as the midpoint of the range between the lower and upper error limits. 
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Any value within the range is equally likely, leading to an infinite number of possible values that can 
be used to calculate spectra. Thus an infinite number of candidate spectral density functions exist. 

Existing spectral analysis methods do not explicitly account for this range of possible 
5 solutions. These methods can be interpreted as using least squared errors or other minimum norm 
methods that attempt to minimize differences between predicted values and measured values. The 
structure of these problems makes them sensitive to such errors when they are solved using 
minimum norm methods, and in particular when solved using least squared errors methods and 
produce an ill-posed problem. These ill posed problems are well known and regularization 

10 techniques are used to calculate reasonable estimated solutions. Regularization methods, including 
Tikhonov regularization (sometimes called ridge regression or Phillips regularization) are used with 
least squared errors fits and augment the least squared errors model with an auxiliary function that is 
extremized (such as minimized). Regularization methods produce stable results, but trade off quality 
of the fitted solution (comparison of the predicted values to the measured data) with some measure 

15 of the regularization function. Making this tradeoff requires using a weighting factor, which is also 
called the regularization parameter, that trades off the weight given to fitting the data with the weight 
given to the regularization function. 

Existing methods for solving ill-posed problems have several deficiencies. They require 
calculating the weighting factor, which is difficult and can have a significant effect on the results. 

20 They fail to explicitly recognize that a solution is only valid when all of the estimated values are 

within the digitization enror bins of the measured values. This deficiency manifests itself in solutions 
that have estimated values that are outside of the digitization bin for one or more measured values, 
and, thus, are not correct. Another deficiency is that in practice the methods rely on using 
minimizing least squared errors, which tends to force estimated values toward the center of the 

25 digitization bins, thus not providing a monotonically convergent means for calculating viable solutions 
by adjusting the regularization parameter. 

Therefore, significant deficiencies exist in the prior art regarding the spectral analysis of 
time-varying signals and signals when it is desired to obtain fine frequency resolution. 

The present invention overcomes problems associated with determining spectral content 

30 from signals where the fundamental frequency of the signal being processed corresponds to a time 
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period longer than what has been measured and/or the underlying spectrum is, at least over a range 
of interest, continuous or substantially continuous (occasionally hereinafter collectively referenced as 
"substantially continuous"). The present invention further enables spectral analysis where the signal 
is nonstationary and where the signal has digitization errors. 
5 More specifically, the present invention is directed to a method and apparatus ("utility") 

involving a time-to-frequency domain transform to define a substantially continuous spectrum or 
other spectrum including non-zero values at irregularly spaced frequency intervals. The invention 
thus allows for spectral estimates for more than a small number of discrete frequencies thereby 
enabling a more accurate representation of certain real spectra and reducing or eliminating spectral 

10 broadening. The invention also allows for examination of a greater range of frequencies of interest 
and provides frequency resolutions that are not limited to a fundamental frequency. In this regard, 
the invention allows for determination of a spectrum based on analysis of a portion of the time- 
based signal that is shorter than a full period of a component of the time-based signal. Moreover, 
the present invention allows for enhanced analysis of received signals that are not stationary, i.e., 

15 that vary with an independent variable such as time. Therefore, the invention has particular 

advantages with respect to medical contexts where the received signal has a changing spectral 
content as a result of interaction of an interrogating signal with moving physiological material such as 
blood flowing through an artery. The invention also enables reduction in digitization or quantization 
errors associated with analog-to-digital conversion. 

20 According to one aspect of the present invention, a utility is provided for processing an 

input signal, such as a medical diagnostic signal, to obtain a spectrum that includes non-zero values 
at irregularly spaced intervals. An associated process involves receiving time-based information 
corresponding to one or more predetermined time intervals of a time-based signal, performing a 
transform on the time-based information to obtain a frequency spectrum including non-zero values 

25 at irregularly spaced intervals, and operating a processor in a signal processing environment for 

using the transform to provide an output based on the time-based signal. The spectrum includes at 
least one pair of successive points (i.e. non-zero values) that has a different frequency spacing than 
at least one other set of successive points. Preferably, the spectrum includes a set of non-zero 
values including a first nonzero value at a first frequency and a second non-zero value at a second 
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frequency greater than the first frequency where the second frequency is a non-integer multiple of 
each frequency of the set of frequencies other than the second frequency. 

According to another aspect of the present invention, a utility is provided for processing an 
input signal, such as a medical diagnostic signal, so as to generate a spectrum based on information 

5 corresponding to a short time interval of the input signal. A corresponding method includes the steps 
of receiving time-based information corresponding to one or more defined time intervals of the time- 
based signal, where the time-based signal includes a component having a period that is longer than 
the time interval, performing a transform on the time-based information to obtain a frequency 
spectrum for the time-based signal, and operating a signal processor in a signal processing 

10 environment for using the transform to provide an output based on a time-based signal. This process 
can be used to obtain a frequency spectrum where the received time-based information 
corresponds to a time interval that is less than half of the period of the signal component. 

According to a further aspect of the present invention, a utility is provided for processing an 
input signal to obtain a substantially continuous spectrum. An associated method includes the steps 

15 of receiving time-based information corresponding to one or more defined time intervals of a time- 
based signal, performing a transform on the time-based information to obtain a frequency spectrum 
for the time-based signal where the spectrum defines a substantially continuous function across a 
frequency range that has non-zero values for a majority of frequencies in the range, and operating a 
processor in a signal processing environment for using the transform to provide an output based on 

20 the time-based signal. Such a substantially continuous spectrum can more accurately characterize 
the real spectrum under analysis, thereby avoiding spectral broadening and other shortcomings of 
certain prior art processes. 

According to a still further aspect of the present invention, a utility is provided for 
processing an input signal so as to account for a digitization or quantization error associated with 

25 analog to digital conversion. An associated method includes the steps of receiving time-based 

information corresponding to one or more defined time intervals of a time-based signal, where the 
time-based signal is an analog signal and the time-based information is digital information, 
performing a transform on the time-based information to obtain a frequency spectrum for the time- 
based signal, where the step of performing a transform involves accounting for a digitization error 

30 associated with a difference between the analog time-based signal and the digital time-based 
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information, and operating a processor in a signal processing environment for using the transform to 
provide an output based on the time-based signal. The invention thus allows for accurate definition 
of a spectrum based on the input signal for certain cases where accurately defining such a spectrum 
may be difficult or impossible due to digitization error. 
5 In accordance with another aspect of the present invention, a utility is provided for use in 

connection with converting a signal from an analog form to a digital form. The utility involves 
receiving an analog signal, defining a number of digitization ranges or bins, and assigning a specific 
value within a bin to a sample related to the analog signal based on a mathematical model where the 
specific value can vary relative to the boundaries of the bin. Preferably, the mathematical model 

10 determines the specific value by solving a constrained optimization problem involving an objective 
function having one or more constraints. For example, such constraints may relate to one or more 
limit values of a bin, a nonnegativity constraint, and a peak count within the bin. For example, the 
objective function may be constrained to be convex and positive within the bin range with slopes 
determined by predefined parameters at or near the bin limits. The constraints may be implemented 

15 as penalty or barrier functions, for example, Heaviside functions. In one implementation such a 
function approximates a Heaviside function and is tapered at an area corresponding to a constraint 
value (e.g., a bin limit value) such that the function is free from singularities (undefined function 
values) at that area. 

In one implementation, the invention is used to determine dimension-related information 
20 regarding a flow channel of a physiological fluid based on analysis of a signal modulated based on 
interaction with the fluid of the flow stream. A corresponding system includes a detector for 
receiving a signal from the fluid and providing an analog electronic signal based on the received 
signal, an analog-to-digital converter for providing a digital signal based on the analog electronic 
signal, a processing means (processor) for first processing the digital signal to obtain spectrum 
25 information corresponding to a frequency spectrum including non-zero values that are not at 
regularly spaced frequencies, and a second processing of the spectrum information to obtain 
dimension-related information for the flow channel. The dimension-related information may be, for 
example, a channel dimension or a volumetric flow rate of the fluid. 
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For a more complete understanding of the present invention and further advantages thereof, 
reference is now made to the following Detailed Description taken in conjunction with the figures as 
described below: 

FIG. 1 Example Spectral Density Function 

FIG. 2 Illustrates an example Gaussian pulse 

FIG. 3 Example Amplitude Function Made from 10 Piecewise Continuous 

Functions 

FIG. 4 Example Amplitude Function Composed of Multiple Functions with 

Overlapping Domains 
FIG. 5 Example Phase Angle with Doppler 

FIG. 6 Example Phase Angle with Doppler when Emitted Frequency Selected to 

Produce Approximately Linear Phase Angle Response 
FIG. 7 Example Phase Angle with Doppler when Emitted Frequency Selected to 

Produce Always Positive and Approximately Linear Phase Angle Response 
FIG. 8 Example sin(Phase Angle) with Doppler when Emitted Frequency Selected 

to Produce Positive Phase Angle Response Close to Zero Phase Angle 
FIG. 9 Example cos(Phase Angle) with Doppler when Emitted Frequency Selected 

to Produce Positive Phase Angle Response Close to Zero Phase Angle 
FIG. 10 Example Phase Angle with Doppler Using 2 MHz Emitted Signal 
FIG. 1 1 Example cos(phase angle) with Doppler Using 2 MHz Emitted Signal 
FIG. 12 Example sin(phase angle) with Doppler with 2 MHz Emitted Signal 
FIG. 13 Six Parameter Distribution for Velocity 

FIG. 14 Quadratic Function Segment Expression for f(t) without Continuity 

FIG. 15 Parameters for Example 

FIG. 16 Example ||A|| Matrix 

FIG. 17 Example ||P|| Matrix 

FIG. 18 Example Envelope Function 

FIG. 19 Spectral Density Component Functions: Used to Calculate Example f(t) 
values 
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FIG. 20 Example Signal and Locations of Measurement Points 

FIG. 21 Example Signal after Adjust for Gaussian Envelope and Locations of 

Measurement Points 
FIG. 22 Example Calculated Results for G(f) and Actual G(f) 
5 HG. 23 Example Calculated Results for F(f) and Actual F(f) 

FIG. 24 Equations for piecewise Continuous Equally Spaced Quadratic Function 
Segments 

FIG. 25 Equations for Piecewise Continuous Equally Space Linear Function 
Segments 

10 FIG. 26 Three Frequency Function with Four Digitization Bins 

FIG. 27 Three Frequency Function with Four Digitization Bins and Least Squared 
Errors Estimate 

HG. 28 Three Frequency Function with Four Digitization Bins and Least Absolute 
Value Errors Estimate 

15 FIG. 29 Three Frequency Function with Four Digitization Bins Least Squared Errors 

and Least Absolute Value Errors Estimates 
FIG. 30 Actual and Estimated G(f) Calculated without Accounting for Digitization 
HG. 3 1 Actual and Estimated F(f) Calculated without Accounting for Digitization 
HG. 32 Bounded Area Measure when Using Piecewise Continuous Quadratic 
20 Function Segments 

FIG. 33 Arc Length Measures for Spectral Density Component Functions Based on 

Piecewise Continuous Quadratic Function Segments 
HG. 34 Bounded Area Measures for Spectral Density Component Functions Based 
on Piecewise Continuous Quadratic Function Segments 
25 FIG. 35 Quadratic Curvature Measures for Spectral Density Component Functions 

Based on Piecewise Continuous Quadratic Function Segments 
HG. 36 Objective Function Combining Measure of Arc Length and Measure of 

Area for both G(f) and F(f) 
FIG. 37 Double Sided Constraint Violation Value Function 
30 FIG. 38 Nonnegativity Constraint Value Function 1 (NCV1) 
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FIG. 39 Nonnegativity Constraint Value Function 2 (NCV2) 
FIG. 40 Tapered Heaviside Gate Function Examples 

FIG. 4 1 Tapered Heaviside Gate Function for Extreme Value of G(f) Less than 
Zero 

FIG. 42 Example G(f) Calculated Using NCV1 Nonnegativity Constraint Value 
Function 

FIG. 43 Peak Count Constraint Value (PCV) Function 
FIG. 44 Basic Calculation Sequence (Page 1 of 3) 
FIG. 45 Basic Calculation Sequence (Page 2 of 3) 
FIG. 46 Basic Calculation Sequence (Page 3 of 3) 
FIG. 47 Actual and Calculated G(f) from Digitized Data (M=50) 
FIG. 48 Actual and Calculated F(f) from Digitized Data (M=50) 
FIG. 49 Actual and Calculated Spectral Density Function from Digitized Data 
(M=50) 

FIG. 50 . Actual and Calculated G(f) from Digitized Data (M=200) 
FIG. 51 Actual and Calculated F(f) from Digitized Data (M=200) 
FIG. 52 Actual and Calculated Spectral Density Function from Digitized Data 
(M=200) 

FIG. 53 Actual and Calculated Cumulative Spectral Density Function from Digitized 
Data (M=200) 

FIG, 54 Block diagram of one embodiment of a system in accordance with the 
present invention 

DETAILED DESCRIPTION 

The interrogating signal emitted from an energy emitter, such as a transducer, can take on 
many forms and EQ. 1 is a general representation of an interrogating signal with frequency/. 
f(t) = E(t)sm(2nft) (i) 

E(t) is a function that describes how the amplitude varies with time. For example, if the 
emitted signal is a pulse then E(t) can be a function that includes a phase shift term that will 
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periodically repeat with period T. An example of a signal that has such a form is a Gaussian pulse. 
The amplitude of a Gaussian pulse will vary with time as defined by EQ. 2 



A Gaussian pulse will vary with time as defined by EQ. 3 and a single Gaussian pulse is 
illustrated in FIG. 2. 



These equations can be used with either pulsed or continuous wave systems. In the context 
of ultrasound medical analysis, the pulsed wave has the advantage of allowing selection of the tissue 
depth to interrogate, which can reduce errors caused by Doppler shifts caused by blood deeper 
than the region of interest, such as the ascending aorta. The pulsed wave approach leads to pulses 
that have a transient response, and this needs to be dealt with, possibly by waiting for the transient 
to pass out of range before collecting data. Another approach could be to model the transient 
explicitly as a pure tone that has changing amplitude or as a pure tone that has multiple amplitudes 
that phase shift with time. Using a continuous wave avoids the transient response issues and leads 
to a pure tone being used to interrogate the tissue. The continuous wave approach, however, has 
the confounding factor of receiving Doppler-shifted signals from regions other than the selected 
depth of the ascending aorta, or whatever other region is of particular interest. Depth dependent 
signals could be inferred from signals that are created from different emitted frequencies (these 
would need to be much different, such as 1 MHz and 3 MHz) such that attenuation effects are 
substantially different. The signal strength differences could be matched against spectra to start 
getting depth-related data. Higher frequencies attenuate more rapidly than lower ones so using 
higher frequencies reduces the interference from deeper tissues, but at the expense of reducing the 
signal from the region of interest. 

As described earlier, the interrogating signal can be used to energize tissue. The resulting 
backscattered signal can be sensed and analyzed to determine useful properties of the tissue that 



E(0 = 5 





(3) 
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backscattered the signal. Such a signal can comprise, for example, the received signal that echos 
back from a multitude of reflectors or backscatterers moving at one or more velocities that have 
been insonified by an ultrasonic interrogating signal. Among the purposes for which the present 
invention can be used is determining the SD functions of backscattered signals from interrogating 
signals that have either constant or time-varying amplitudes, including time varying signals that 
include one or more pulses, repeating sets of pulses, or any other amplitude function that varies 
with an independent variable. As discussed above, time is an example of such an independent 
variable. 

For illustrative purposes, the following description and derivations will use as an example 
signals such as those that may occur when ultrasonic signals are used to interrogate tissues in a living 
organism. This illustration is not limiting in that the present invention does not require using 
interrogating signals or a response from such a signal. The measured data used by the present 
invention may be acquired using any suitable means that eventually produces two or more digitized 
data, including transducers followed by subsequent processing or data input or read from one or 
more sources of data that may have been previously stored, such as tables or graphs stored as hard 
copy or using electronic, magnetic, or optical means. 

Let A(tf) be the SD function (in this case the function is the same as an amplitude density 
function). The general form of the signal for which the SD function is to be obtained is EQ. 4. 



As described in more detail below, the present invention determines SD functions for 
processes with two or more independent variables, such as signals that vary with time, g(t), when 
the SD function depends on at least one independent variable,/, that has one or more nonzero 
values over one or more parts of the interval \f L ,f H ], where f L mdf H are the lower and upper limits 
off. EQ. 4 defines an amplitude function, A(t,f), that varies based on the values of two 
independent variables, t and / For this example, the variable t represents time and the variable/ 
represents frequency. The Aft f) function can be separated into functions of /and of L 




A(/,/)-B(OBCO (5) 
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When EQ. 5 is substituted into EQ. 4 the result is a general expression for the time varying 
signal g(t) for which the SD function B(f) is sought. 



g(0 = | E(r)B(flsin(27c/r + df (6) 



EQ. 6 can be rearranged to become EQ. 7 

g(0= E(OB(/)(sin(2 7t/Ocos( ft/)) + cos(2tc/0 sin(<K/))) ( 7 ) 

Hie following terms are defined for coefficients. 

F(/)=B(/)sin(<K/)) (8) 
G(/) = B(/)cos(<|>(/)) (9) 

Substituting these equations for the coefficients into EQ. 7 leads to EQ. 10 

g(0 = E(0 F(/)cos(2rc/0 + G(/)sin(2 7t/0<*/ (10) 
\ 

Distributing the integral over the sum of the sine and cosine terms leads to EQ. 11. 



g(0 = E(0 



j ?(f) cos(2 nf t) df+ j Q(f) sin(2 nft)df 



(11) 



Define spectral density component functions (SDC functions) as functions from which the 
SD function can be calculated. F(f) and G(f) are SDC functions. As is well known, F(f)andG(f) 
15 can be used to calculate the spectral density function B(f) using EQ. 12. 
B(/) = VF(/) 2 +G(/) 2 (12) 

The amplitude functions can also be used to calculate the phase shift for each frequency 
using EQ. 13. 

<K/)=arctan(^j (i 3 ) 

20 SD functions and SDC functions both characterize how the amplitude of a function varies 

with respect to an independent variable and are collectively referred to as amplitude functions. For 
example, F(f) characterizes how the amplitude of F(f) cos(2 Tift) varies with frequency. Either the 
SD function itself or SDC functions are to be estimated using the means described below. 
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The SD function or the SDC functions, such as F (f) and G(f), may be represented using 
various functions. The functions used to represent amplitude functions will generally have 
parameters that, when selected correctly, will lead to useful estimated values that are calculated for 
g(t). Amplitude functions may be represented by any functional form for which parameters can be 
5 determined that produce estimated values for g(t) that differ from the actual values of g(t) in a 
manner that meets predetermined objectives. An example of such a predetermined objective is 
having the errors between the measured values and estimated values of g(t) being small enough that 
the estimated values for g(t) provide useful approximations of the actual values of g( t). Useful 
values for g(t) occur when the parameters specified for the functions that represent the SDC 

10 functions can be used to calculate the SD function of g(t). 

One aspect of the present invention involves structure and methodology for using values of 
g(t) and E(t) y either based on measurements, simulations, or other suitable means, to determine the 
specific values of the parameters that characterize the SD function or one or more SDC functions 
when such amplitude functions are composed of one or more substantially continuous functions or 

15 when such amplitude functions use values for one or more independent variables, such as frequency, 
that lead to non-orthogonal basis functions constituting g(t). Using continuous functions to represent 
one or more amplitude functions over at least part of the domain of an independent variable (such as 
time) is particularly advantageous when a multitude of signals with various frequencies combine to 
produce a composite signal that is the measured signal g(t). Such a composite signal occurs, for 

20 example, when an interrogating signal backscatters from the multitude of particles traveling over a 
wide range of velocities as they are being carried in a flowing fluid. An example is the signal 
produced when an interrogating ultrasonic signal backscatters from the blood cells traveling through 
a living organism's aorta. 

In contrast to the present invention, Fourier methods represent the amplitude functions using 

25 discrete values of frequencies that are selected as integer multiples of a fundamental frequency. The 
frequencies used in Fourier methods are carefully selected such that the resulting basis functions are 
both orthogonal and normal (orthonormal). As described earlier, using discrete frequencies, 
particularly when they are restricted to being harmonics of a fundamental frequency, to approximate 
continuous functions leads to significant inaccuracies, such as spectral broadening. The present 

30 invention improves upon Fourier methods by using frequencies that can be selected substantially 
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without limitation as to which frequencies are selected and also by allowing the use of continuous 

functions to represent continuous amplitude functions. 

In the present invention, one or more functions are selected to represent the SD or SDC 

functions. The function used to represent an SDC, such as F(f)> may be zero. The functions used 

5 to represent the amplitude functions may be discrete, such as one or more vectors of discrete 

frequencies starting at any frequency of interest and spaced at any, possibly varying and nonuniform, 

frequencies from each other (in contrast to the evenly spaced harmonic frequencies using in Fourier 

methods). The functions representing amplitude functions in the present invention can be orthogonal 

or non-orthogonal. The functions representing amplitude functions can be either normal or not 
i 

10 normal. The function or functions used to represent amplitude functions may be one or more 

continuous, not continuous, or piecewise continuous functions. The functions representing amplitude 
functions may have one or more intersecting, including overlapping, domains. An overlapping 
domain between two or more functions is an interval where the domains of two or more functions 
intersect at more than one value. FIG. 3 illustrates an example where multiple piecewise continuous 

15 functions that intersect at only their endpoints represent an amplitude function. FIG. 4 illustrates an 
example where multiple functions with overlapping domains model an amplitude function. One 
aspect of the present invention is that functions representing the amplitude functions may use one or 
more of the following functions: splines, B -splines, polynomials of order 1 or greater, trigonometric 
functions, trigonometric polynomials, exponential functions, hyperbolic functions, Heaviside 

20 functions, functions that approximate Heaviside functions (discussed in detail below), or any another 
functions that can be used to construct one or more bases that span one or more predetermined 
vector spaces that can approximate or define a discrete, non-continuous, continuous, or piecewise 
continuous function. For example, the piecewise continuous segments illustrated in FIG. 3 are 
composed of quadratic polynomials and Heaviside functions. As another example, FIG. 4 

25 illustrates how B-splines can combine to represent an amplitude function. 

FIG. 4 also illustrates another aspect of the present invention that provides for 
computational flexibility. The independent and dependent variables can be offset (as illustrated in 
FIG. 4 where the frequency axis is offset by 2 MHz) or scaled to facilitate calculating the 
parameters of the functions that represent amplitude functions. 
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Using piecewise continuous functions or functions with overlapping domains is of particular 
interest because parameters for function segments defined over discrete intervals can be efficiently 
determined using the present invention. The union of the function segments can represent 
complicated amplitude functions. FIG. 3 is a graphical example of a piecewise continuous function. 
5 The points where the individual function segments of the piecewise continuous function intersect are 
knot points and these are illustrated by the hollow dots in FIG. 3. 

Another aspect of the present invention involves structure and methodology employing a 
calculation module that uses one or more objective measures to determine the values of the 
parameters in the functions that represent the amplitude functions or SD function. One example of 

10 such an objective measure is using goodness of fit. Other examples of objective measures are 

minimizing the length of the amplitude functions, SD function, or SDC functions, minimizing the 
squared area or absolute value of the area under the amplitude functions, SD function, or SDC 
functions, and minimizing one or more measures of the slopes or changes in slopes of the amplitude 
functions, SD function, or SDC functions. 

15 One example of an objective measure using goodness of fit is minimizing, or approximately 

minimizing, the sum of the squared errors between measured values of g(t) and those predicted 
using amplitude functions. Another example of an objective measure of goodness of fit is 
minimizing, or approximately minimizing, the sum of squared errors when the values of one or more 
of the parameters in one or more of the amplitude functions or the SD function are constrained to 

20 have values with predetermined properties. Such predetermined properties may include that the 
parameters are allowed to only have particular values, such as being non-negative, or that the 
parameters have a particular relationship, such as the mathematical combination of two or more of 
the parameters (such as one or more of their sum, difference, product, or quotient,) always be 
positive or that the mathematical combination always falls within a certain range, such as being 

25 approximately equal to zero. Another example of an objective measure of goodness of fit is 

minimizing, or approximately minimizing, the sum of the absolute values of errors, including the 
possibility that the values of one or more of the parameters in one or more of the amplitude 
functions or the SD function are constrained to have values with predetermined properties. As 
before, such properties include having one or more parameters being restricted to particular values, 

30 such as being non-negative, or that the parameters have a particular relationship, such as the sum of 
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two or more parameters always be non-negative or that they approximately equal zero. These 
aspects of the invention are covered in more detail below. 

Referring to EQ. 11, in accordance with a further aspect of the present invention, optimized 
results can beobtained by setting the values for f L and/„ respectively low enough and high enough 
that they cover the full range of frequencies that constitute the spectrum of g(t). The values for/ L 
and/tf may be determined concurrently with determining the values of the parameters for the 
functions representing the SDC functions, such as G(f) and F(f), or they may be selected prior to 
determining the parameters for the functions representing G(f) and F(f). Selecting/ L and f H ahead 
of time is often possible. For example, the conditions that produce g(t) are generally known. A 
predetermined lowest frequency may be based on the frequency of an interrogating signal interacting 
with one or more moving media. For example, the lowest frequency can be predetermined based 
on the frequency of the interrogating signal and the lowest frequency expected from the moving 
media after including factors such as measurement errors or noise that can affect the frequencies 
produced by the interaction of the interrogating signals with the moving media. 

Another aspect of the invention involves using single continuous functions that span the 
entire domain of frequencies as amplitude functions or as SD functions. An example of such a 
function is EQ. 14. This equation characterizes the amplitude functions that occur during laminar 
fluid flow with measurement noise. This equation would become B(f) in EQ. 8 and EQ. 9. 

(14) 



(c + a)f e 



I 1 + * 



-2- 
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(/•+/ 8 ) 2 C*+/ S )2a 2 ; 
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-dx 



The emitted frequency sent out as the interrogating signal is f e and the tissues move at 
various velocities to produce the reflected frequencies,/. The amount of tissue at each velocity 
produces the relative amplitudes. The parameters are: 

1. F the fraction of the insonified region that has flowing liquid. 

2. t the shape exponent for the flow region 
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3. cs s : the standard deviation for the velocity in the slow/nonflowing region's reflectors. 

4. f e : the frequency of the interrogating signal. 

5. a: the maximum value for velocity that is achieved at the time that a measurement is 

made, 

6. \i s : the mean velocity of the slow moving/nonflowing region's reflectors. 

7. or ; the standard deviation of the measurement error associated with the flowing fluid. 

8. c: the speed of sound. 

A further aspect of the present invention involves using such a single continuous amplitude 
function in combination with a function that represents the phase shift, <p(f), see, e.g., EQ. 8 and 
EQ. 9, to produce single continuous functions that represent amplitude functions that span the entire 
domain of frequencies. An additional aspect of the present invention involves using such a single 
continuous amplitude function, e.g., in combination with functions that represent the sin((p(f)) and 
cos( $f)) in EQ. 8 and EQ. 9 to produce single SDC functions that span the entire domain of 
frequencies. A further aspect of the present invention relates to representing <p(f), sin((p(f)), or 
cos( <p(f)) using either constants or polynomials with degree of 1 or greater. A further aspect of the 
present invention involves selecting the frequency of the interrogating signal used with Doppler 
measurements of velocities such that <p(f), sin( $/)), or cos(<p(f)) may be represented using either 
constants or a linear function. 

A further aspect of the present invention involves characterizing the phase shift when using 
Doppler measurements of fluid flows using EQ. 15 to represent phase angle. 



(J) = -arctan 

wh 
ere 



(15) 

\ 



2/71 cosj^— j -fn+c bf sin\^- jcosjj-^— j-/7te J ' 



(J) = phase angle in radians 
f = frequency in Hz 
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D = depth of region from which Doppler-shifted signals are being generated 
c = speed of signal, such as speed of sound 
x = attenuation exponent 
b = attenuation constant 

For the case where x=l EQ. 15 simplifies to become EQ. 16 



<|) = arctan 



-2 n sin 


(2' D7t } + C *co(2 /D ^ 


- e be 


2 71 cos 


' 2 fDn> 


+ c b sin 




- 2 % e 



(16) 



The attenuation constant, b, can be calculated from the attenuation factor using EQ. 17 
b = .5000000000 10 " 7 (3 ln( 10 ) (17) 
where P is the attenuation factor in dB. 

FIG. 5 illustrates how the phase angle can change with frequency. As is known, selecting 
the emitted frequency determines the frequencies of the returned signal. As previously mentioned, 
one aspect of the present invention involves selecting the emitted frequency such that the returned 
signals are confined to an approximately linear region of the phase angle response. For example, in 
FIG. 5 selecting an emitted frequency that leads to return frequencies around 2.079 MHz would be 
in a region where phase angle varies close to linearly with frequency, whereas selecting an emitted 
frequency that leads to return frequencies about 2.17 MHz would be in a region where phase angle 
does not vary linearly with frequency. More specifically, if 2.078 MHz is selected as the emitted 
frequency and blood flow in the ascending aorta of an adult human is being measured using an 
interrogating signal beamed down the axis of the aorta, then frequencies in the range of 
approximately 2.078 to 2.093 MHz will be returned. FIG. 6 illustrates an example phase angle 
response over this frequency range. In this example, selecting 2.078 MHz leads to phase angle 
responses that may be modeled as a linear function. 

FIG. 5 illustrates that half of the phase angle response will be over frequencies that lead to 
positive phase angles. Selecting an emitted frequency such that the phase angle response is always 
positive is beneficial. When the phase angles are always positive the amplitude response function 
F(f) of EQ. 1 will always be positive. If the phase angle becomes negative the sin((p(f)) term will 
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be negative, leading to negative values over at least part of the range of F(f). Having only positive 
phase angles avoids this condition. Calculating correct values for F(f) and G(f) is made easier 
when it is known that F(f) and G(f) are always positive. Selecting the emitted frequency for 
Doppler velocity measurements such that the phase angle response is always positive is included in 
5 accordance with another aspect of the present invention. 

Another aspect of the present invention involves selecting the emitted frequency so that the 
phase angle response is always positive and in a region where the phase response is approximately 
linear. FIG. 7 illustrates an example in accordance with this aspect of the present invention. 

A further aspect of the present invention involves selection of the emitted frequency so that 
10 at least a portion of the phase angle response is close to zero. 

FIG. 8 and FIG. 9 illustrate how thus selecting the emitted frequency yields sin( <p(f)) and 
cos( <p(f)) functions that can be approximated. sin( <p(f)) can be approximated in many cases using 
a linear function. The range of cos( (p(f)) is small. Depending upon the precision needed it may be 
approximated as a constant (such as 0.98 in this example), as a linear function, or as a higher order 
15 polynomial. A linear approximation would pick up most of the variation and be adequate for many 
applications. 

Compare the results illustrated in FIGS. 7, 8 and 9 with what happens if the emitted 
frequency is not selected according to this aspect of the present invention. FIGS. 10, 1 1, and 12 
illustrate the results from using an emitted frequency of 2 MHz for the conditions of this example. 
20 The linear response and small domains that occurred when the emitted frequency was selected in 
accordance with this aspect of the present patent are not present in FIGS. 10, 11, and 12. 

When the frequency distribution of a spectral distribution has been produced in a manner 
that relates frequency to a calculated variable then the distribution of the calculated variable can be 

25 determined. For example, if the frequency distribution has been obtained using Doppler then the 
velocity at any frequency can be calculated using EQ. 18. The frequency data or data for a 
calculated variable can be used to calculate the values of one or more parameters that characterize 
a distribution. Such distributions are particularly useful when the methods of the present invention 
are used to calculate continuous distributions. Such a distribution could be, for example a Gaussian 

30 distribution that is characterized by its mean and standard deviation. EQ. 19 is an example of a 
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more complicated distribution that has six parameters. FIG, 18 is an example of such a distribution 
when the calculated parameter is velocity. 



c (-/+/,) 




Other aspects of the present invention involve structure and methodology for determining 
the parameters of an SD function, or the SDC functions, from signals that have been distorted by 
digitization errors. When signals are measured and converted from analog to digital form the analog 
values are converted to numeric (digital) values with a limited precision. Such limited precision 
leads to numeric values that approximate the actual values of the signal and differences between the 
actual values and the approximate values are errors caused by the digitization process. In 
particular, the present invention includes structure and methodology for calculating SD functions, or 
the functions that can be used to determine the SD functions, when using digitized data by 
constraining the calculated functions so that they have predetermined properties. Such 
predetermined properties can include having at least one, and typically the majority or 
approximately all, of the values of g(t) estimated using the calculated functions falling within the 
ranges determined by the digitization process. Such predetermined properties can also include 
requiring SD functions, or the functions that can be used to determine the SD functions, having 
predetermined properties related to their geometries or of the parameters that define the functions. 
For example, such methodology includes having functions that are the minimum measures of slopes 
or rates of changes in slopes, the lengths, or squared areas enclosed by the functions. Another 
example is having functions that maximize, or minimize, the magnitudes of one or more of the 
parameters that characterize the functions. Another example is having functions that minimize 
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measures of length and squared areas enclosed by functions and that maximize the magnitudes of 
one or more of the parameters that characterize the functions. These aspects of the invention are 
described more completely later. 

To illustrate the aspect of the present invention that characterizes amplitude functions using 
multiple functions, consider the following description of a general process for estimating amplitude 
functions. Amplitude functions, such as those of EQ. 8 and EQ. 9, can be approximated arbitrarily 
closely as the sum of other functions. The following expressions define the approximations of G(f) 
and F(f) as the sums of N piecewise functions. 



The N piecewise functions defined over each of N intervals will collectively span the domain 
of /that is of interest. N functions exist for G(f) and N functions exist for F(f), and they may or may 
not have overlapping domains. Each of the Gfj) and F/f) will be defined over a domain of /with 
lower and upper limits. For the case of using piecewise continuous functions the union of the 
domains will be continuous and the individual domains will not overlap except at the points of 
intersection. 

Substituting EQs. 19 and 20 into EQ. 1 1 and distributing cosine and sine terms over the 
sums leads to EQ. 21. 



(19) 



N 



(20) 



(21) 
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Changing the integral of the sum to a sum of integrals leads to EQ. 22. 
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g(0 = E(0 



F.(f) cos(2 nft)df 



(22) 
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,(f)sm(2nft)df 
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Each of the piecewise functions is defined over part of the domain off. Indexing the limits 
of integration for each of the integrals in the sums leads to EQ. 23. 



g(0 = E(0 



N r J 



(23) 



F.(f) cos(2% ft)df 



S G/f) $m(2% ft)df 



10 



£(f j is the envelope of the interrogating signal. If the signal being analyzed has not been 
modulated by an interrogating function then setting the envelope function E(t)=l will decompose the 
received signal without removing the effects of the interrogating signal. If E(t) is not constant then 
methodology in accordance with the present invention can separate the effects of Eft) while 
decomposing the received signal into, for example, SD or SDC functions. 

g(t) is the received signal and measured values for it are known. Dividing both sides of EQ. 
23 by Eft) produces EQ. 24. 



s(Q 
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F,(f) cos(2 Kft) df 
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(24) 

\ 



£ G(/)sin(2 7c/l)d/ 



Measured values from g(t) are the data used to calculate F(f) and G(f). 
15 Use the following substitutions in EQ. 24.: 

g(0 



E(0 =f(0 (25) 
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^Mj- H j (26) 



(27) 



and obtain EQ. 28. 



f(/) = 



f H. 
N r J 

If F i 

JmlJ ■ 



\ / fi, 

N n i 



F,(f) cas(2nft)df 



E f W)si 
7-1 J, 



(28) 
sin(2 7t/r)d/ 



10 



As stated earlier, the piecewise continuous functions may be polynomials, such as linear or 
quadratic functions. Each of the function segments will have its coefficients, such as the slope and 
intercept for linear functions. For illustration of the principles of the present invention, the following 
example uses quadratic functions and derives expressions that may be used to solve for the 
coefficients that specify each of the quadratic function segments. 

Let the following expressions specify the piecewise quadratic function segments: 

G.(f) = a j f + b j f+c j (29 ) 

F(f)=x J f+y J f+z J (30) 

Substituting these expressions into EQ. 28 leads to EQ. 31. 



f(0 = 



S (*/ + cos(27t ft) df + X {a.f + b J f+c J )sm{2nft)df 



(31) 

\ 



15 This equation is linear in the coefficients x p y j7 z jt a p b p and c Jt and the integrals have 

analytical solutions. FIG. 14 has the analytical solution. This equation is a general expression. No 
restrictions have been placed on the interval limits [Lj, Hj[ so they do not need to be continuous and 
they could overlap. Therefore, the function segments could be discontinuous. No constraint has 
been placed on the slopes of function segments matching if function segments are continuous. The 

20 current formulation, without constraints, is broadly applicable. 
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The unknowns to be solved for are the coefficients x j9 y p z P a p b p and c } for j=i . . N. 
Represent these unknowns as the vector |d|. Once values have been estimated for |d| then 
estimated values for f(t) can be calculated using EQ. 31. The values of |d| can be solved for using 
one or more objective criteria, such as the estimated values for f[t) having the least sum of squared 
5 errors compared to the measured values for /ft). EQ. 31 is linear in x p y p z p a p b p and cj so 
minimum norm methods, such as least squared errors, may be employed to solve for the 
unknowns. The number of frequency intervals and the lower and upper frequencies for each 
function segment are values that are selected, and thus, are known. Similarly the times when the 
measurements are taken are known. These values are inserted into the equation in FIG. 14 and the 
10 result is a set of coefficients that multiply each of the x j% y p Zp a p bp and c } . This set of coefficients 
forms a rectangular matrix ||A||. Designate the vector of measured values of f(t) as |f|. The resulting 
matrix equation is EQ. 32. 

||A|||d| = |f| (32) 

Matrix solution methods known to those skilled in the art can be used to solve for |d|, which 
15 contains the values for x p y p z p a p b p and c } . One such method is using the Moore-Penrose 

generalized inverse of ||A||. Designating ||P|| as the Moore-Penrose generalized inverse of ||A[| leads 
to EQ. 33. Singular value decomposition is another method that can be used. 
IPII|f| = |d| 03) 

For the situation where the measurement times and frequency intervals are predetermined 
20 and reused multiple times only one matrix inverse needs to be calculated, which offers considerable 
economy of calculation. The matrix ||A|| is based on the predetermined [L p Hj\ and predetermined 
measurement times and is calculated in advance. The Moore-Penrose generalized inverse, |[P||, is 
also calculated in advance and stored for later use, such as in the nonvolatile memory of a 
computation means such as using a hard disk or other magnetic media or on a CD ROM or other 
25 optical storage media, or nonvolatile RAM or other electronic storage media. After j|A|| and ||P|| 
have been calculated and stored the values for ||P|| can be retrieved and used to calculate ]d| 
whenever such calculations are desired for a vector of measured values of f(t), |f| using EQ. 33, 
which leads immediately to values for the coefficients xp y Jy zp a p bp and cj. 

As described below, the formulation becomes substantially simpler when the lower and 
30 upper bounds for each interval, [Lp Hj\ all have the same width, but they do not need to have the 
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same width or be equally spaced. No restrictions are placed on the measurement times being 
equally spaced. The flexibility available for selecting the bounds for each interval and the 
measurement times allows selecting nonuniform intervals, such as using Chebyshev points or 
expanded Chebyshev points. 

Methods other than the least sum of squared errors can be used to calculate the values of 
the unknown coefficients . For example, the least sum of absolute value of errors (LAV) or 
nonnegative least squares (NNLS) can be used. As is known to those skilled in the art, LAV can 
be implemented such that constraints are placed on the values of the variables, as is also the case 
with NNLS. Such constraints are particularly useful when linear function segments are used 
because these readily lend themselves to constraints that prevent one or more of SDC functions, 
such as G(f), from having any negative values. 

For the case where spectra are continuous, such as the Doppler spectra from flowing 
blood, the jth function segment needs to connect to function segment j-L Requiring continuity 
means that the endpoint of the {j-l)th function segment must be the starting point for the jth function 
segment, which leads to the lower bound of the jth segment being equal to the upper bound of the 
(j-l)th segment. Therefore, EQ. 34 holds for piecewise continuous functions. 

L r H j-i (34) 

As demonstrated below, this restriction means that the value for c. j s determined by the 

values of a j {9 b $ _ x% c.^and Therefore, ^ is not an independent variable when the 
quadratic function segments are piecewise continuous. 

Spectra from processes such as flowing blood are smooth. Smooth transitions between 
quadratic function segments requires that the slope at the end of the (j-1 )th segment match the slope 
of the jth segment. As will be shown below, this constraint removes the set of bj from being 
independent variables. 

The following discussion uses as examples the 0 ■> b., c. t and the same discussion also 

applies to the *j> z. , Include an offset for the frequency at the beginning of the interval over 
which a function is being used. The equation for the jth function segment that goes through a 
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specific point is derived as follows. Designate the specific point for the lowest frequency and the 
amplitude at the lowest frequency by using the ordered pair ( Zy P^ ). Pj is the beginning amplitude 
for the;ffc function segment. The jth function segment has its own coefficients a v b- p and intercept 

5 The starting point for the )th function segment has the same coordinates as the endpoint for 

the (j-l)th function segment Designate the endpoint for the (j-l)th function segment as 

( 1? e Jmm x y is the ending amplitude for the (j-l)th function segment. The starting point for 

the jth function segment is the endpoint of the j-lth function segment. The following equation gives 
the amplitudes for the jth function segment of the SDC function when using piecewise continuous 
10 quadratic functions to represent the SDC function. 

Gj-ajV-H^f + btf-Hj^ + Cj (35) 

However, c.-e, _ so 
j j-i 

G r a J (f-H j _ l )\b J (f-H j _ { )^e J _ l (36) 

15 For the special case when j=l the starting point has an amplitude of 0 at the frequency f lb , 

where f lb is defined earlier. Clearly, a recurrence exists that gives the e k . The basic equation is: 

e r a j {H r H j _ l )\b j {H r H j ^)^e j ^ (37) 

20 This equation gives the logical result that the height of the endpoint of the jth function 

segment is the sum of the heights of all of the preceding endpoints plus the change in height that 
occurs during frequency interval j. Using the preceding equations to eliminate Cj from EQ, 35 leads 
to EQ. 38. 
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(38) 

\j>ml ^ J 

These equations that tie the endpoints together establish continuity between the piecewise 

functions. The piecewise functions are now piecewise continuous. The values for the H } are 
known because they are selected a priori. In this set of equations the only free variables are the 
first N-l ^ and the N a.. 

Constrain the slopes of the j-lth and the jth segments to match where the segments joint. 
This constraint leads to the starting slope of the jth function segment being equal to the parameter 

b. The bj for each segment depends only upon the values for the coefficients for the previous 
function segments. The end result is that each bj depends only on the previous values of the a.. 
The equation for & is 



6=2 
j 



( J \ 
Y(H -H J a , 



(39) 



For the first function segment the initial slope is zero because it starts tangent to the f-axis. 
Therefore, ^ = 0. This situation becomes clear by examining the case where j=l. 
G r a x (J-H Q f (40) 

Using EQ. 39 to eliminate b from EQ. 38 leads to EQ. 41 . This equation is the general 
expression for the jth function segment for piecewise continuous quadratic function segments. 



The endpoint of the final function segment is at frequency H . The SD function needs to be 
zero at this point and the slope also needs to be zero. These requirements also require the SDC 
functions, such as G(f), to have zero amplitude and slope at H N . Two boundary conditions exist 
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at H N , but only one free variable, a N _ { , remains, therefore need to remove another degree of 
freedom. Set the equations for the boundary conditions that equate the slope and amplitude at the 
end of the Nth segment equal to each other (they both equal 0) and it becomes clear that a , is 
not a free variable, either. Therefore, when have N function segments only N-2 values of |a| need to 
be found. 

The equations become considerably simpler is equally spaced function segments are used. 
Assume that each of the function segments covers a frequency interval of equal width. Select a 

value for the lowest basis frequency, H Q , and then index by up from it for other frequencies in 
the basis. 

With equally spaced frequency intervals have 

H q-r H g -2= H & (41) 

And the general expression for the jth quadratic function segment becomes EQ. 42. 

(42) 



J 5 



X (-(2/>-l)ff 5 + 2/-2if 0 )o M + fl (f-ff o -0'-i)/f 8 ) 
\p =1 ) 



The unknown variables are the a } . Represent these variables as the vector |a|. When 
impose continuity and matching slopes at the knotpoints where the function segments join the values 
for the (N-l)th and Nth values of |a| stop being independent variables. The expressions for these 
variables are 

N-2 

L 

P 



a N - x =lL%(P-m (43) 



V-2 \ 

X a (p-N+l)\ (44) 

j» i J 



/N-2 



The equations that give a N _ x and a N confirm that the last two coefficients are completely 
specified by the first N-2 coefficients and, thus, are not independent variables. These values lead to 
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the last two function segments working together to make the final amplitude of G(f) become zero 
and with zero slope with smooth transitions at the knotpoints of the final two function segments. 
Expressions for the (N-l)th and Nth function segments in terms of the first N-2 elements of 

(45) 

N-2 

I(ff 8 (-(2;-l)ff 5 + 2/-2ff 0 )fl.+ (-H 0 -(W-2) H & +/) 2 a } {j - N) ) 

|a| are 



(46) 

N-2 

G N = I (-«■(-/+ NH 5 + H 0 ) 2 U-N+l)) 
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Similar relationships exist for F(f) y Fj and |d|, including 



(47) 

2 



F ff)=H h ^H2 r l)ff 8 + 2/-2^l + (-// fl -(/'-l)// 8+ /) x 

J 



(48) 



N- 1 
N-2 



£ {H s {-(2j-l)H s +2f-2H 0 )x+(-H 0 -{N-2)H & + rfx.U-N)) 



10 



T/-2 



15 



/W=S(^W^+^) CH*+D) (49) 

Substituting the preceding equations into EQ. 22 and integrating and substituting to remove 
the N-lth and Nth aj and x } leads to EQ. 50. This equation is the final equation for the function that 
approximates g(t) when using equally spaced quadratic function segments that are piecewise 
continuous and constrained so that the slopes at the knotpoints are equal and that slopes and values 



20 



(50) 



at the beginning and end of the SD function are zero. EQ. 50 is a general formula for calculating 
g(t) that is linear in the variables x. and a jm 

The times when measurements are made are known, as are the starting basis frequency, H Q 
and the width of each frequency interval, H* After these substitutions are made the result is a set of 

equations consisting of coefficients on the unknown N-2 unknown variables a. and the N-2 

unknown variables *. . The coefficients can be extracted and put in a rectangular matrix ||A||. ||A|| 

will have M rows and 2(N-2) columns where: 

M = number of measurements, with the measurements being t h i = 1 . . M 

N = number of function segments (which is the same as the number of basis frequencies). 
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Set ||A|| so that the first N-2 columns contain the coefficients for aj and the last N-2 
columns contain the coefficients for the x Jt The equation for the first N-2 columns of j|A|| is EQ. 51. 
The equation for the last N-2 columns of ||A|| is EQ. 52. These equations assume that the envelope 
function is known, which is often the case, such as if the interrogating function is a continuous sine or 

(51) 

A w =^E(l,)(2co8(2n(H 0 + yH 6 )i / )-2w)s(2n(H 0 +y- l)» a )f,) 

+ 2tN-j-l)cos(2nt 9 {H 0 + NH 6 )) 

+ 2 (N-j) cos(2 n t. (H 0 + (N- 2) ff,)) 

+ {-4N+4j+2)c<x(2nt l (H Q +(N- 1) ff 5 ))) /(n 3 1') 

(52) 

^, i + ^-2 = i E ^/) (" 2 J" 1 ) sin(2 K f. (f/ 0 -I- 
-2sin(2 7t(tf 0 +jff 5 )g^ 
+ (4iV-4j-2)sin(2 7cr.(tf 0 +(iV- 1)#*)) 
+ 2 sin(2 7C (ff Q + U~ 1 ) H 6 ) f,)) /(n 3 1) 

has other known characteristics. If the signal being analyzed has not been modulated by an 
interrogating function then set the envelope function E(t)=l and the method will decompose the 
complete signal and not attempt to remove the effects of the interrogating signal. 

Designate the values of the (currently unknown) N-2 values of a. as |a| and designate the 

(currently unknown) N-2 values of *. as |x|. Designate the vector |d| as being |a| stacked over |xj, 
thus |d| contains the unknown variables a p j = 1 . . N-2 followed by the unknown variables 
j = . . N-2. The values |dj can solved as described earlier, using, for example, the Moore-Penrose 
generalized inverse. 

FIG. 15 through FIG. 23 illustrate the method just described using an example. The 
example uses a small number of measurements, 20, and seven quadratic function segments. FIG. 
15 lists the parameters used to calculate the ||A|| matrix shown in FIG. 16. The Moore-Penrose 
generalized inverse of ||A|| is shown in FIG. 17. Only ||P|| is needed to calculate the values of |d| 
from measured values of f(t). The 20 values of g(t) used in the example were calculated using the 
G(f) and F(f) shown in FIG. 19. A successful calculation will produce estimated G(f) and F(f) close 
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